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Abstract 

We present the LWL formula which represents the long wavelengh limit of the solutions 
of evolution equations of cosmological perturbations in terms of the exactly homogeneous 
solutions in the most general case where multiple scalar fields and multiple perfect fluids 
coexist. We find the conserved quantity which has origin in the adiabatic decaying mode, 
and by regarding this quantity as the source term we determine the correction term which 
corrects the discrepancy between the exactly homogeneous perturbations and the k — > 
limit of the evolutions of cosmological perturbations. This LWL formula is useful for 
investigating the evolutions of cosmological perturbations in the early stage of our universe 
such as reheating after inflation and the curvaton decay in the curvaton scenario. When 
we extract the long wavelength limits of evolutions of cosmological perturbations from the 
exactly homogeneos perturbations by the LWL formula, it is more convenient to describe 
the corresponding exactly homogeneous system with not the cosmological time but the 
scale factor as the evolution parameter. By applying the LWL formula to the reheating 
model and the curvaton model with multiple scalar fields and multiple radiation fluids, we 
obtain the S formula representing the final amplitude of the Bardeen parameter in terms 
of the initial adiabatic and isocurvature perturbations 
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1 Introduction and summary 



Recently we come to be required to investigate the evolution of cosmological perturbations 
in the very early universe [26]] , [6]] . According to the inflationary scenarios and the curvaton 
scenario, in the early universe the wavelength of cosmological perturbations responsible for 
the present cosmic structures such as galaxies and clusters of galaxies is much larger than 
the horizon scales. Therefore the methods for researching the cosmological perturbations 
on superhorizon scales have been sought. In this context, Nambu and Taruya pointed out 
that there exists an LWL formula representing the k — > limit of the cosmological pertur- 
bations in terms of the exactly homogeneous perturbations [32]]. Soon later in the multiple 
scalar fields system [15]] [30]], the complete LWL formulae were constructed. Since the 
evolution equations of the corresponding exactly homogeneous universe look simpler than 
the evolution equations of cosmological perturbations, the LWL formula brought about 
great simplification. In addition, the viewpoint that the evolutions of the cosmological 
perturbations on superhorizon scales are governed by the stability and instability of the 
corresponding exactly homogeneous universe [9]], [10]] is useful for physical interpretation. 
In this context, the phase space of the corresponding exactly homogeneous system was 
investigated in detail and the role of the fixed points in the phase space in the stability 
and instability of cosmological perturbations was discussed [10] ] . For these reasons, the 
LWL formula was used for investigating the evolution of cosmological perturbation on su- 
perhorizon scales by several authors [H], [30]], [32]], [9]], [TO]]. In this paper, in section 
2, in order to investigate the evolutionary behaviors of cosmological perturbations during 
reheating and the curvaton decays, we construct the complete LWL formulae for the most 
general system where the multiple scalar fields and the multiple perfect fluids coexist. 

In the early universe, the cosmological perturbations on superhorizon scales responsible 
for the cosmic structures and the CMB temperature anisotropics experience the reheating 
and/or the curvaton decays. In these processes, the multiple scalar fields such as inflatons 
and curvatons oscillate coherently, gradually decaying into radiation fluid. By replacing the 
oscillating scalar fields with the dust fluids, the evolution of the cosmological perturbations 
during reheating after the inflation [8]] and in the curvaton decay [20]] were investigated. 
These authors treated the system dominated by dust-like scalar field fluid and radiation and 
investigated the influence of the entropy perturbation originating from the multicomponent 
property to the evolution of the total curvature perturbation variables such as the Bardeen 
parameter. The purpose of these analyses was to determine the initial perturbation of 
the present Friedmann universe in terms of the early stage seed perturbation. Although 
it was shown partially that this replacement is physically reasonable [H], [8]], we are 
required to treat the decaying oscillatory scalar fields directly. In fact, the instabilities 
characteristic to the rapidly oscillating scalar fields were pointed out [3J], [5]], [37]], [TU1 ]. 
[361 ] and investigated [33J ] , [TO] ] , [3T] ] , [UJ ] . In order to treat the oscillatory scalar fields 
directly, the action angle variables were introduced and the averaging method representing 
the averaging over the fast changing angle variables was applied [9]], [TO]]. In this paper, 
in order to investigate the evolutionary behaviors of cosmological perturbations during 
reheating and/or the curvaton decays, in section 3 the action angle variables and the 
action angle perturbation variables are introduced, and in section 5, the averaging method 
was applied to the decaying oscillatory scalar fields. 

In the papers [9]], [TO]], the averaging method was used to investigate the corresponding 
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exactly homogeneous system, and by the LWL formula the evolution of the cosmological 
perturbation in the long wavelength limit was constructed from the corresponding exactly 
homogeneous perturbation. By using the LWL formula and by applying the averaging 
method, in single oscillatory scalar field system [14 ] and in nonresonant multiple oscillatory 
scalar fields system [9] ] it was shown that the Bardeen parameter is conserved, and in 
resonant multiple oscillatory scalar fields system [10]] it was shown that the cosmological 
perturbation including the Bardeen parameter can grow. In this paper, in section 6, 7, by 
using the LWL formula and by applying the averaging method to the decaying oscillatory 
scalar fields, we construct S formulae representing the final amplitude of the Bardeen 
parameter in terms of the initial seed adiabatic and entropic perturbations, in the reheating 
and in the curvaton decays, respectively. 

The organization and the summary of the paper is explained as follows. In section 
2, we construct the LWL formula as for such scalar-fluid composite system based on the 
the philosophy of the paper [15] ] . The discrepancy exists between the evolution equations 
of the cosmological perturbations in the k — > limit and the evolution equations of the 
exactly homogeneous perturbations because the former contains fc 2 $ = 0(1) terms and 
the latter does not, therefore this discrepancy should be corrected by the correction term 
which contributes the adiabatic decaying mode, but any general methods for determining 
such correction term have not been presented yet, and only in the multiple scalar fields 
system such correction term was determined. We show in the /c — *> limit the existence 
of the conserved quantity which has origin in the adiabatic decaying mode and which is 
related with /c 2 $. By regarding this conserved quantity as the source term, and obtaining 
the special solution A^, we correct the exactly homogeneous perturbation A" and we obtain 
the complete LWL formula A = + in the most general scalar-fluid composite system. 
In section 3, we point out that it is more appropriate to use the scale factor a rather than 
the cosmological time t as the evolution parameter when we use the LWL formula. As for 
the scalar quantity T, we use the perturbation variable DT, D is the operator which maps 
the exactly homogeneous scalar quantity T to the gauge invariant perturbation variable 
representing the T fluctuation in the flat slice. D defined in this way can be interpreted as a 
kind of derivative operator. In fact, the exactly homogeneous part (DT)* can be expressed 
as the derivative of T with respect to the solution constant with the scale factor a fixed. 
In order to investigate the exactly homogeneous system containing oscillatory scalar fields, 
we use the action angle variables I a , 8 a . By using D defined in this way, we can define the 
action angle perturbation variables DI a , D9 a whose exactly homogeneous parts are given 
as the derivatives of I a , 9 a with solution constant C with the scale factor a fixed. When 
we use the derivative operator D and the LWL formulae, it is essential to use the scale 
factor a as the evolution parameter. In section 4, we apply the LWL formulae to the non- 
interacting multicomponents system and discuss the long wavelength limit of the evolution 
of the Bardeen parameter. In section 5, we apply the averaging method by which the 
system is averaged over the fast changing angle variables to the decaying scalar fields which 
have been discussed in the reheating model and the curvaton model. By evaluating the 
corrections produced by the averaging process and the errors produced by the truncation of 
the sufficient reduced angle variables dependent part, the validity of the averaging method 
is established. In section 6, 7, we apply the LWL formula and the averaging method to 
the interacting multicomponents model such as the reheating model, the curvaton model, 
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respectively. We assume that the multiple scalar fields and the multiple radiation fluid 
components exist. In these models, we construct the S formulae representing the final 
amplitude of the Bardeen parameter in terms of the initial adiabatic and isocurvature 
perturbations. In our previous paper [10 ], the evolutionary behaviors of cosmological 
perturbations in the early universe where multiple oscillatory scalar fields interact with 
each other have been investigated. This S formula give the information about how the 
cosmological perturbations which grew in such early era are transmitted into the radiation 
energy density perturbations through the energy transfer from the scalar fields into the 
radiation fluids. We present the necessary condition for the initial entropic perturbations 
produced in the early era to survive until the late radiation dominant universe. Section 
8 is devoted to discussions containing non-linear generalization of our LWL formailism 
and comment of the case where the decay rate depends on other physical quantities. In 
appendices, the proofs of the propositions presented in section 5 and the evaluations of the 
useful mathematical formulae used in section 6 are contained. 

In this paper, we consider the case where the homogeneous scalar fields obey the phe- 
nomenological evolution equations as 

dU 

+ 3#0+ — + S = 0. 1.1 
ocp 

The interactions between scalar fields are described by the interaction potential U, while 
the interaction between scalar fields and fluids are described by S. This analysis includes 
the well known case S = Tip [11] ], [21] ] and the general case where S is an arbitrary 
analytic function of 0, which was discussed in the paper [35] ] but whose perturbations 
have not been investigated yet. Another supplemental purpose of our paper is to present 
the evolution equations of cosmological perturbations corresponding to the homogeneous 
system with various source term S especially dependent on <fi. 

The notation used in this paper are based on the the review [12]] and the paper [TBI]. 



§2 Derivation of the LWL formula 

We give the definitions and the evolution equations as for the background and the pertur- 
bation variables. Based on these notations, in the most general model where the multiple 
scalar fields and the multiple perfect fluid components interact, we give the LWL formula 
representing the evolutions of the perturbations variables in terms of the exactly homoge- 
neous solutions. 

We consider perturbations on a spatially flat Robertson- Walker universe given by 

ds 2 = -(1 + 2AY)dt 2 - 2aBY j dtdx j 

+a 2 [(l + 2H L Y)5 jk + 2H T Y jk ]dx i dx k , (2.1) 

where Y, Yj and Yj k are harmonic scalar, vector and tensor for a scalar perturbation with 
wave vector k on flat three-space: 

Y Y 3 :=^Y, Y jk := (±5 jk - ^Y (2.2) 
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By using the gauge dependent variables 1Z and a g representing the spatial curvature per- 
turbation and the shear, respectively: 

K:=H L + ^H T , a g := ^H T - B, (2.3) 
we can define two independent gauge invariant variables: 

A:=A-(^)\ *:=K (2.4) 

In order to define the matter perturbation variables, we will consider the scalar quantity 
perturbation variables generally. As for covariant scalar quantity T = T + STY, we define 
the gauge invariant perturbation variable representing the T fluctuation in the flat slice: 

DT := 5T - ^-TZ, (2.5) 
H 

Next we consider the covariant scalar quantity T 2 whose background quantity is the time 
derivative of T: T. The extension of T into the covariant scalar quantity T 2 is not unique. 
For example, 

/ \ I - ~-|l/2 

T 2 = sgn (doTj \ -r u V,TV v T\ , (2.6) 

and 

T 2 = n^T, (2.7) 
where is an arbitrary vector field satisfying 

n% = -1, (2.8) 

have the same T as the background part. But these defferent T 2 's give the unique pertur- 
bation part: DT 2 = (DT)' — TA. Therefore we can define DT by 

DT := (DT)' — TA. (2.9) 

We consider the the universe where the scalar fields a , (1 < a < N s ) and the fluids 
p a ,P a (1 < a < Nj) coexist, whose energy momentum tensor is divided into A = (S, f) 
parts where S represents the multiple scalar fields, / represents the multiple fluids. The 
energy momentum tensor of / part are further divided into individual fluids parts a. On 
the other hand, the energy momentum tensor of S part cannot be divided into individual 
scalar fields parts a, since the interaction potential U contains the terms consisting of plural 
scalar fields a : 

fM = (fc) s +(f{j) = (f^ s + J2T^, (2.10) 

a 
a 

where the energy-momentum transfer vector Qa\i is defined by 

V.T^ = Qa/j, = Qau^ + /a m , (2.12) 
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where is the four velocity of the whole matter system and the momentum transfer f Afl 
satisfies u^f Alx — 0. For the scalar perturbation, the energy- momentum tensor and the 
energy-momentum transfer vector of each individual component are expressed as 

f% = -(pa + 6 Pa Y), (2.13) 

T% = a(p A + P A )(v A -B)Y j , (2.14) 

TL = (P A 6i + 6P A Y6i + n TA Yj), (2.15) 

and 

Qao = -[Qa + (QaA + 5Q a )Y], (2.16) 

Qaj = a[Q A (v-B) + F cA }Yj, (2.17) 

where p A , P A and Q A are the background quantities of the energy density the pressure 
and the energy transfer of the individual component A, respectively. The anisotropic 
pressure perturbation Ut A and the momentum transfer perturbation F cA are already gauge 
invariant. As for the scalar quantities T = (p A , P A , Q A ), we use DT as the gauge invariant 
perturbation variables. As for the gauge invariant velocity perturbation variable, we use 

Z A :=n-^-(v A -B). (2.18) 
The energy-momentum tensor of scalar fields part is given by 

(f£) s = • - l - (V A • V A + 2C/) 62. (2.19) 
Since divergence of the energy momentum tensor is given by 

(V M T£) s = (nj> a - SCj V,0 a , (2.20) 

in order that the phenomenological equations of motion of the scalar fields become 

dU 

□0a - -y- = S a , (2.21) 

0<pa 

we assume that 

•V,V,x V (2.22) 



Qu) 



s 

By using the scalar fields background variables (f> a , <p a , S a and the corresponding pertur- 
bation variables D(p a , D<j) a , DS a , the background part of the fluid variables are given by 

Ps = \($ l + U, (2.23) 

Ps = \($f~U, (2.24) 

h s = (0) 2 , (2.25) 
Qs = S-j>, (2.26) 
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and the perturbation part of fluid variables are given by 



(Dp) s = 


d Ps . D <P + dPs ■ D<P, 
d<P V d<j> 


(2.27) 


(DP)s = 


o<p d(f) 


(2.28) 


(hZ) s = 


-H<j> ■ D(j), 


(2.29) 


(n r ) s = 


o, 


(2.30) 


(DQ)s = 


-S-Dj>-j>- DS, 


(2.31) 


(aF c ) s = 


S a (-kD(f) a - jj<fiaZj . 


(2.32) 



When the source of the scalar field (j) a , S a is given as functions of the covariant scalar 
quantities T and T 2 whose background part is T, that is S a = S a (T,T 2 ), DS a is given by 

DS a = H ■ DT + ^ ■ DT. (2.33) 

In such case, (DQ)$ can be written as 

(DQ)s = 8 -§f-DT+ d -^-Dt, (2.34) 

which is assumed from now on. In the same way as the individual components T^ u , as 
for the total energy-momentum tensor = Avi we can define the gauge invariant 
perturbation variables such as Dp, DP, hZ and Hp. From (12. 101) . (12.111) . we obtain the 
background equations as 



p 


= Ps + ^pa, 

a 


(2.35) 


p 


= p s +J2 p ^ 

a 


(2.36) 


h 


= hs + ^K, 

a 


(2.37) 





= Qs + ^ Q a , 


(2.38) 
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and perturbation equations as 



Dp 


= D Ps + ^Dp a , 

a 


(2.39) 


DP 


= DP s + J2 DP «i 

a 


(2.40) 


hZ 


= {hZ) s + ^h a Z a , 

a 


(2.41) 


n T 


= (IT T ) 5 + ^IT TQ , 

a 


(2.42) 





= {DQ) s + J2DQ a , 

a 


(2.43) 







(2.44) 



This Z is known as the Bardeen parameter [2]] [12]], [23]]. In the long wavelength limit, the 
Bardeen parameter is conserved in the case where the entropy perturbations are negligible. 
But in various systems it was reported that the entropy perturbations cannot be neglected 
0], [51], P], [10], so in the present paper we will investigate the evolutionary behavior of the 
Bardeen parameter more carefully. Until now, as for the gauge invariant scalar quantity 
perturbation variables, we use D. But traditionaly most scalar quantity perturbation 
variables have been written without using D: 

Y a := £>0 a , p a A ga := Dp a , P a U La := DP a , Q a E ga := DQ a . (2.45) 

This Y a has been called the Sasaki-Mukhanov variable [28] ] [22] ] . 

In terms of the gauge independent variables defined above, we give the evolution equa- 
tions of cosmological perturbations. From (12.211) . the background and the perturbation 
parts can be written as 

dU 

+ 3#0 + — + S = 0, (2.46) 

L 1 {DT,A) = ~D<l>-*£*, (2.47) 
er a z H 



where 



d 2 U ■ ■ dU 

L^DT, A) = {D<py + ZH{D<py + -tttk-D^ + DS-(pA+ 2(— + S)A. (2.48) 



As for the fluid components, V^T^ = Q au gives the background equations as 

p a = -3Hh a + Q a , (2.49) 

and the perturbation equations as 

k 2 

L 2a (DT, A) = ~^jjh a ($ - Z a ) , (2.50) 

/,mZ " ^ + 3h a Z a + h a A + DP a - lu Ta = -jF ca + %-Z, (2.51) 



HI 3 k H 
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where 

L 2a (DT, A) = (D Pa y + 3HDp a + 3HDP a - Q a A - DQ a . (2.52) 
G% = k, 2 T^ gives the background equations as 

H 2 = yp, (2.53) 

p = -3Hh, (2.54) 

H = --(l + w)H 2 , (2.55) 



and the perturbation equations as 



L 3 (DT,A) = 2p— (2.56) 



3a 2 H 2 

K 2 ^ k 2 

— n T — 

3 a 

3 



L 4 (DT, A) = -— n T - , (2.57) 



^+-(1+«;)Z = 0, (2.58) 

• 4 + j(l*)-F a2n - < 2 ' 59 » 



where 



L 3 (DT, A) = 2pA + Dp, (2.60) 
L 4 (DT,A) = HA + 2HA - — (Dp + DP). (2.61) 

The dynamical perturbation variables are classified into two groups, that is, what has 
analogy with the exactly homogeneous perturbations and what is not related with the 
exactlty homogeneous perturbations at all. The dynamical perturbation variables of the 
former type are DT representing the scalar quantity T = (p, P, <p, Q, S) perturbation in 
the flat slice, DT and the metric perturbation variable A. The dynamical perturbation 
variables of the latter type are the Newtonian gravitational potential $ and Za, F c a, ^-ta 
which have vector or tensor origin. In the above L^i — 1, • • •, 4) equations, the former 
type dynamical perturbation variables are contained in the left hand side while the latter 
type perturbation variables are collected in the right hand side. The exactly homogeneous 
perturbations DT^ and A^ corresponding to DT and A, respectively are constructed as 

{DT)> ■= (w),-^' (262) 

A> := (2.63) 

*' - !g); < 264 » 

where C is the solution constant of the background solution and the subscript t implies 
that the derivative with respect to C is performed with the cosmological time t fixed. On 
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the other hand, the dynamical perturbation variables of the latter type such as $, Z4, 
F c a and Uta do not have exactly homogeneous counterparts. The evolution equations 
of cosmological perturbations containing Li(i = 1, • • -,4) have analogy in the exactly 
homogeneous perturbation equations. In fact, the variations of the exactly homogeneous 
equations (12^61) . (|H9j) . fl233|) and ([235]) give 

Li(DT*,A*) = (i = l,...,4), (2.65) 

respectively. The only difference between the exactly homogeneous perturbation Li(i = 
1, • • •, 4) equations and the actual k 7^ cosmological perturbation Li(i = 1, • • •, 4) equations 
is that fc 2 $ terms exist in the latter but k 2 Q terms do not exist in the former. Then the 
effect of the source term fc 2 $ is corrected in the following way. In performimg the correction 
process, it is important to notice that the source terms fc 2 $ can be represented in terms of 
conserved quantity which has origin in the universal adiabatic decaying mode. In fact, as 
for / defined by 

/ 1 \ k 2 
f = a 3 H [A + - A g )= — a$, (2.66) 



3H 



using (12301) . fl236|) . (12371) yields 



= -a 3 H 2 wU T + ^ak 2 (l + w)Z. (2.67) 

When we assume that for k — ■> limit 

IT T -»• 0, kZ -> 0, (2.68) 

are satisfied, the quantity / is conserved, whose value is written as c. Therefore for k — > 
limit, 

fc 2 $^^c = 0(l). (2.69) 

This expression of $ is well known as that of the universal adiabatic decaying mode [TBI]. 
In the Li(i = 1, ■ ■ -,4) equations containing DT, A, the Newtonian potential $ appears 
only in the form fc 2 $, that is, accompanied by k 2 . When we assume that DT = 0(1), A = 
0(1), /c 2 $ behaves as 0(1). Since in the linear perturbation, the scale of the perturbation 
variables is arbitrary, the fact that $ = 0(l/fc 2 ) does not imply the breakdown of the 
linear perturbation. If one want to get $ = 0(1), one simply assumes that DT = 0(k 2 ), 
A = 0(k 2 ). But as explained later, we cannot assume that $ is vanishing, since c defined 
by (12.691) must satisfy the constraint (12.801) . Therefore in the k ^ limit where (I2.68p . 
(1239|) are satisfied, (l2~4Tj) . (12301) . d236|) . (12371) can be written as 



L la (DT,A) = 


30a 
3 C ' 


(2.70) 


L 2a (DT,A) = 


3h a 

3 c ' 
a 6 


(2.71) 


L 3 (DT,A) = 


2p 


(2.72) 


L 4 (DT,A) = 


3H 

3 c - 

a 6 


(2.73) 
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It can be verified that above four sets of equations (I2.70p . (I2.7ip . (I2.72p and (12.731) are 
satisfied by 



A 



DT 



3 -a 

2 V 
f 

H 9 > 



\ , 9 



where 



dt- 



(2.74) 
(2.75) 

(2.76) 



to 



This special solution for A = (DT,A) is written as A . Since the variation of the exactly 
homogeneous solution A 9 satisfies (EMJ, the general solutions of (12770]) . (12771]) . (T2772]) . 
(I2.73P A = (DT, A) can be expressed as 



A = A 6 + A\ 



(2.77) 



The perturbation equations except Lj(i = 1, • • •, 4) equations have vector origin, that 
is, they are derived from the space component of the Einstein equations. Therefore these 
perturbation equations do not have any analogy with the exactly homogeneous perturba- 
tion equations. As explained in the paper [151 ]. these perturbation equations determine the 
evolutions of the dynamical perturbation variables which have vector or tensor origin, that 
is, which have no correspondence with the exactly homogeneous pertubations, or give the 
constraint which should be satisfied in order that the exactly homogeneous perturbations 
become the k — > limit of evolutions of cosmological perturbations. Therefore (I2.5ip . 
(I2.58P can be interpreted as the decision of the evolution of the variables Z a which is not 
related to the exactly homogeneous solution at all in terms of A = (DT, A), the constraint 
to the exactly homogeneous perturbations, respectively. Integrating (I2.5ip yields 



H 



C a + I dta 3 ( -h a A 

to 



DP a — —F ca + Q^-Z, 
k H 



(2.78) 



By summing (12.781) with respect to all the fluid components, we obtain 



(hZ) ] 



H 



H 



V] C a + f dta 3 (-h f A - (DP) f - S 



2 a 3 
3i7 



pA+a 3 



2 a 3 
3H 



pA + a 3 <p ■ Dcf) 



(2.79) 



Therefore (I2.58P gives the constraint between C a , c defined by (I2.69p . and 2N$+Nf solution 
constants of the exactly homogeneous perturbation as 



0. 



Integrating (I2.59j) gives 



H 

a 



C t - dtaA^j 



(2.80) 



(2.81) 
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where the first term containing C t is well known universal adiabatic decaying mode [15]] 
and by comparing with (I2.69P we obtain 

C t = ^c. (2.82) 

If we assume c = 0, since (12.801) gives one constraint relation, we obtain 2N$ + 2Nf — 1 
solutions and the Newtonian potential is obtained by (I2.8ip with C t = 0: 

H f 

$ = / dtaA, (2.83) 

a J t 

If we assume that c is nonvanishing, since c = 0(1), A 0(1), therefore 

^0(k% (2.84) 

we obtain the k — > limit of the universal adiabatic decaying mode [TBI]: 

-fc 2 $ -> —c, (2.85) 
3 a 

which is consistent with (I2.69p . Then we have obtained the long wavelength limit of all 
the solutions to the evolution equations of cosmological perturbations. 



§3 Use of the scale factor as the evolution parameter 

As the gauge invariant variable representing the fluctuation of the scalar quantity T, we 
adopt DT defined by (12.51) which represents the T fluctuation in the flat slice, since it is the 
easiest to see the correspondence with the exactly homogeneous perturbation of T. While 
until now we described the exactly homogeneous variables as functions of t, C where t is the 
cosmological time and C's are solution constants, we can describe the exactly homogeneous 
variables as functions of a, C where a is the scale factor. For an arbitrary scalar quantity 
such as p, P, S, Q, <fi, from (I2.62p . (I2.64p . (DT)^ can be written as the partial derivative of 
the corresponding exactly homogeneous scalar quantity T with respect to solution constant 
C with the scale factor a fixed: 

c*T =(§).- («) 

Since 

mi-- mi-- 

this property of D also holds as for the time derivative of the scalar quantity T: 

d 



gc T) ■ (3.3) 



a 



Therefore the operator D defined by (12.51) can be interpreted as a kind of derivative oper- 
ator, that is, the derivative with respect to the solution constant C with a fixed. Because 
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of this derivative property of D, as for the scalar quantities T = (p, P, S, Q) which are 
functions of 0, <fi, we can understand 

dT dT 

DT = — -D<f ) +- I .D^ (3.4) 

easily. The contribution to DT from the adiabatic decaying mode is given by the same 
form as that of DT: 

f -\ b f 

[DT) =-g, (3.5) 

where g is defined by (I2.76p . 

As seen from the above duscussion, in order to derive the long wavelength limit of 
cosmological perturbations from the corresponding exactly homogeneous system by using 
the LWL formula, it is more appropriate to use as the evolution parameter the scale factor 
a than the cosmological time t. For example, as for the scalar-fluid composite system, 
the corresponding exactly homogeneous expressions are obtained by solving the first order 
differential equations setting (fi a , p a := a 3 (f> a , p a as independent variables and the scale 
factor a as the evolution parameter: 



a 3 

% (3.6) 
replacing H with the right hand side of the Hubble law: 



d A 


1 Pa 


a-T<Pa = 
da 


Ha 3 ' 


d 


a 3 dU 


a— p a = 
da 


Hd<p a 


d 


-3h a + 


a— p a = 
da 



H 2 = - 



2a 6 



(3.7) 



While under use of t as the evolution parameter our system is the constrained system with 
the Hamiltonian constraint (13.71) . under use of a as the evolution parameter our system 
becomes the unconstrained system with respect to independent variables a , p a := a 3 a , 
p a . The corresponding first order perturbation variables are D<f) a , P a := a 3 D<f) a , Dp a . 

For some time, we consider the system consisting of multiple scalar fields cf) a only. Since 
the evolutions of 4> a , p a := a 3 <p a can be described in terms of the Hamilton equations of 
motion, the evolutions of the corresponding perturbation variables Y a = D(f) a P a := a 3 D<j) a 
can also be written in terms of the Hamilton equations of motion: 



dY a OH dP a dH 



dt dP a ' dt dY a 



(3.1 



whose Hamiltonian is given by 



H = -Lp a P a + ^-V ab Y a Y b + ^J b P a Y b , (3.9) 
2a- 5 2 2H 

d 2 u 3k 2 ■ ■ k 2 fdu ■ ■ du\ k\ . 

Vab = + ~ Mb + m[Wa h + ^WJ + ^ 5ab - (3 - 10) 
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When we discuss the quantization of the fluctuations, the other set of the canonical per- 
turbation variables Y a := Y a , P a := a 3 Y a has been used [231 ]• But in the viewpoint of the 
LWL formalism, the set of the canonical variables Y a := D(j) a , P a := a 3 D<fi a is more natural 
than the set of the canonical variables Y a := Y a , P a := a 3 Y a , because the long wavelength 
limit of the former set is generated from the derivative of the homogeneous variables <j) a , 
p a := a 3 <p a with respect to the solution constants with the scale factor fixed. The con- 
nection between the old canonical variables Y a := Y a , P a '■= a 3 Y a and the new canonical 
variables Y a := D(j) a , P a '■= a 3 D(f) a are given by the canonical transformation defined by 
the generating function: 

W = Y a P a + l — j>akYoY h . (3.11) 

4 p 

When we treat the oscillatory scalar fields, the action angle variables I a , 9 a are useful 

mi ma]: 



7,3/2 



2I a a 
COS 6 a , 



Pa 



,3/2 



a/ 2m a I a sin 9 a , 



(3.12) 



where m a is the mass of the scalar field 4> a - The action angle variables obey the evolution 
equation as 

a 3 / 2 



d 

a-rla 
da 

d ft 

da 



a 3 dUu 



+ 



H 36 a H 



11 a 



m 

3/2 



sin 9 a S a + 3/ a cos 29 a , 



ma a 3 dU int cr 

H H dl a H ^J a 



cos 9 a S a - -sin20 a . 



(3.13) 
(3.14) 



In order to investigate the cosmological perturbations in the universe containing oscillatory 
scalar fields, by using D defined in the above we define the action angle perturbation 
variables DI a , D9 a starting from Y a '■— D<f) a , P a '■= a 3 D(j) a . In the LWL formalism, the 
perturbation variables corresponding with the action angle variables I a , 9 a are DI a , D9 a 
defined by the following expressions: 



D 



D 



a 3 / 2 

,3/2 



21, 



a n 

- cos U a 



-a 



2m a I a sin 9 a 



where D in the right hand side is interpreted as 



_d_ 



(3.15) 



(3.16) 



The expressions obtained from variations of I a , 9 a with a fixed in the previous papers [9] 
], [10] ] are the long wavelength limits of DI a , D9 a defined by (13. 15ft . In fact, the LWL 
formulae 



DL 



D9 n 




9, 



(3.17) 
(3.18) 
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where g is defined in (12.761) hold. While the third term in the right hand side of (13.171) 
appears because of the scale factor a dependence of the transformation law from <p a , <p a 
to I a , 9 a , the Jj parts of (I3.17P and (13.181) reflect the fact that D is the derivative operator 
with respect to the solution constant with the scale factor a fixed. 

In order to solve the dynamics of the system containing the oscillatory scalar fields, we 
are required to perform the averaging over the fast changing angle variables 9 a [9]], [TU1 ]. 
If we use the cosmological time t as the evolution parameter, our system is a constrained 
system. Therefore we must check that our averaging procedure is consistent with the 
constraint and this process is rather cumbersome. But if we use the scale factor a as the 
evolution parameter, our system becomes unconstrained system, so the definition of the 
averaging procedure becomes rather simple. 

We can conclude that use of the scale factor a as the evolution parameter brings about 
the two merits. One is that it becomes easier to see the correspondence between the exactly 
homogeneous solution and the long wavelength limit of the first order perturbation and 
that the LWL formulae become more simple. The other is the more simple definition of 
the averaging process. 

We consider the evolution of the Bardeen parameter Z. Following the paper [3H ], 
we define ( as the gauge invariant variable representing the curvature perturbation in the 
uniform density slice: 

C := -— Dp = K - — bp. (3.19) 
P P 

From (12.561) (12.581) we can see that the Baredeen parameter Z and ( are closely connected 

as 

2 1 k 2 
9 1 + w a z H z 

whose k — > limit is 

2 1 r 

Z = C-~— (3-21) 

where c is constant related with the adiabatic decaying mode defined by (I2.69p . While the 
Bardeen parameter Z is expressed as the weighted sum of Z s , Z a which do not related 
with the exactly homogeneous quantity at all and whose evolution is written in the rather 
cumbersome integral form (I2.78p . ( is connected with the corresponding exactly homoge- 
neous quantity and (/ evolution can be written in terms of the derivative of the total energy 
density p with respect to solution constant. Then we consider much easier Q evolution. 

In the paper [191]. as the nonlinear generalization of the Bardeen parameter (, ((t,x) 
was introduced. When P = P(p), ((t,x) is reduced to 

1 M*' 33 ) do 

C(f , x) = \na(t,x) + -j + p ( ) • ( 3 - 22 ) 

In fact, the first order quantity of ((t, x) is given by 

tf,.)=^ + ;M (3-23) 
a 3 p + P(p) 

which agrees with the Bardeen parameter (. In the viewpoint of the LWL formalism, we 
adopt the zero curvature slice da(t, x)/dx l = 0. We assume the equation of state P r = p r /3, 
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since the final state of reheating and of the curvaton decay is radiation dominant. In this 
case, the above ((t, x) is reduced to 



((a,x) = - In p r (a, C(x)). 



(3.24) 



p r (a,C(x)) is the expression of p r obtained by solving the locally homogeneous system 
(see the separate universe approach [34|]) with use of a as the evolution parameter. The 
solution constants C(x) have spatial dependence. By considering C(x) = C + 5C(x) and 
expanding with respect to 5C(x), we can obtain the perturbation of an arbitrary order up 
to the decaying modes of perturbations. For example, the n-th order perturbation is given 
by 



-T 



d" 



dC ai dC a2 ■ ■ ■ dC an L4 



- In p r (a, C) 



5C ai (x)SC a2 (x)---SC an (x). 



(3.25) 



When we obtain the exactly homogeneous exression p r = p r (a,C), we can know the long 
wavelength limits of perturbations of arbitrary orders. Later we determine the exressions 
Pr — Pr{ a i C) in the reheating and in the curvaton decay. 

In the reheating and in the curvaton decay, 5C(x)'s are given by the action angle 
variables in the initial time: 



5I a (x) := 5I a (a = a , x), S9 a (x) := 59 a (a = a , x). 



(3.26) 



We discuss how to determine the statistical properties of 5I a (x), S9 a (x). 5I a (x), 59 a (x) 
are given at the time when the slow rolling phase ends and the coherent oscillation begins. 
At this time, <fi a = 0, and as for the perturbations 



x 



d A ke lk - x e a (k), 5<p a (x) = 0, 



where e a (k) is the Gaussian random variable satisfying 
< e a (k)e b (k') >= P a (k)5 ab 5 3 (k + k') 



k := Ifcl. 



By solving I a 9 a in terms of <p a <p a and by Taylor expanding, we obtain 



SlJx) 



x 



a 3 m a (f) a 5(f) a (x) + -a 3 m a [5(j) a (x)Y 
0, 



(3.27) 

(3.28) 

(3.29) 
(3.30) 



where we use <f> a = 0, 5c/) a (x) = 0. The above fact that 89 a (a = ao,x) = does not 
imply 69 a (a, x) = 0, since 9 a (a) is a function depending on not only 9 a (a ) but also / a (oo)- 
Therefore it is often important to consider the role of the perturbations of angle variables. 

Determining the many point correlation function of fluctuations is reduced to the eval- 
uation of 

< e ai (k 1 )e a2 (k 2 ) ■ ■ ■ e an (k n ) > . (3.31) 
This quantity can be determined by applying the differential operation defined by 



exp 



£ /*fcP.(*>^ 



(*)M-*) 



e=0 



(3.32) 
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§4 Application of the LWL formula to the non-interacting 
multicomponent system 

Based on the results obtained in sections 2, 3, we consider the long wavelength limit of 
the evolutions of cosmological perturbations in the universe consisting of multiple cosmic 
components. For simplicity, we consider the case where each component do not interact 
with each other. 

We consider the w a fluid where w a is constant and which does not interact with other 
components Q a = 0. p a is solved as 

A a 

Pa = a 3(l+w a ) ' ( 4 - X ) 

where A a is a solution constant and by differentiating with respect to A a with the scale 
factor a fixed, we obtain 

<*W = s^y. (4 - 2) 

where SA a is a perturbation solution constant corresponding with A a . Therefore we obtain 

Al = = const. (4.3) 

Next we consider the case that the oscillatory scalar field <fi a does not interact with other 
cosmic components, that is U- lxtk = 0, S a = 0. Since in such case the right hand side of 
(I3.13P is oscillatory function depending on the angle variable 6 a with vanishing mean value, 
by taking the averaging over 9 a , we obtain 

I a = A a , (4.4) 

where A a is constant. The estimate about the effects of the oscillations due to the fast 
changing angle variables 6 a was discussed in the previous papers [9] ] [10] ] . Therefore by 
taking the derivative with respect to solution constant A a with the scale factor a fixed, we 
obtain 

DI\ = 5A a} (4.5) 
where 5A a is a perturbation constant corresponding with A a . Since 

m a I a m a DI a 

Pa = — 5-, Dpa = o — , (4.6) 

a 6 a 6 

where the above second expression is given by the D operation to the above first expression, 
we obtain 

A L = ^ --1T = const. (4.7) 



A n 



As seen from the above expressions, as for the oscillatory scalar field <p a the energy density 
p a and the energy density perturbation p a ^ ga behave like those of the dust fluid. We can 
summarize that for non- interacting cosmic components, Aj^, A$ ga are conserved. 
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As for the multicomponent non-interacting fluids system, evolution is given by 

C=(E;&)/(»&+^)- <«■ 

We can see that (} is exactly conserved for the adiabatic growing mode defined by 

5A a 1 



A a l+w 



a independent. (4.9) 



For some time, we consider the two components system consisting of dust and radiation. 
In this case, (14.81) is reduced to 

* = <"°> 

where the suffix d, r imply dust, radiation, respectively. We obtain 

A A - l JA± (4 11) 

in the limit a — > 0, a — > oo, respectively. We adopt different, more physical paramerization: 

5A d = 3£A d - ^ <L4 r = 4£A r . (4.12) 

£ represents the adiabatic growing mode and rj represents the isocurvature mode defined 
by 

CL = 0, Si = 3 -Al - Al = const =: rj. (4.13) 

Then we obtain 

CL = CL - lv, (4-14) 



which is the famous formula. [ 

Although in the paper [3U], (14.101) has already been derived essentially without using 
the LWL formula, our result is more rigorous in the point that we treat the contribution 
from the adiabatic decaying mode characterized by c defined by (12.69P more appropriately, 
while the paper [34 ] simply assumes that k 2 § is vanishing. 

§5 Application of the averaging method to the decay- 
ing scalar fields 

We derive the evolution equations of the multiple scalar fields decaying into the multiple 
radiation fluids. By solving these evolution equations and taking the exactly homogeneous 
perturbations, we can obtain the information of the evolutionary behaviors of cosmological 
perturbations during reheating and in the curvaton model. We assume that the source S a 
is given by 

Sa = T a a - (5.1) 
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We nondimensionalize the dynamical quantities as 



— -> a, t- -> io, — ^ Pa, — -> t/"int = (5.2) 
«o A) Po Po 



(5.3) 



and the parameters as 

m r 

where z/ is the small parameter implying the ratio of the interaction energy to the free part 
energy defined by 

Po ■= —a-- (5-4) 
Then we obtain the dimensionless parameters 

« pj /2 v^ro , , 

e:=^= , 7 — — 1^-, (5.5) 

which imply the ratio of the Hubble parameter to the mass of the scalar fields H/m, the 
ratio of the decay rate to the Hubble parameter Y/H at the initial time a = ao, respectively. 
By using the above dynamical variables and parameters, the evolution equations can be 
expressed as 

a^-h = ---^ TJ ^^--'yY a ^f-(l-cos29 a ) + 3I a cos2e a , (5.6) 
d 1 m a 1 a 3 dU int 1 T a 3 , . 

«^ = ijn + 7jr^-^jr^ 2 ^-r m26 - (5 - 7) 

a ^°" Q = ^ r aa m a / a (1 - cos 26> a ) , (5.8) 
where a a is defined by 

P« = -J, (5-9) 

a 4 

and T aa is the decay rate from the scalar field I a to the radiation component a a and 
therefore r a is given by 

Y a = J2 r *a- (5.10) 

a 

In this paper, we investigate the evolutionary behavior of cosmological perturbations during 
the period when the decay rates from the scalar fields to the radiation fluids are large 
compared to the interaction between the scalar fields, that is 7 ^> u/e, while in the paper 
[TU1 ] the evolutions of cosmological perturbations during the period when the interaction 
between the scalar fields is dominant u/e ^> 7 were discussed, which is thought to give the 
initial conditions for the present studies. 

Next we show that there exists a transformation such that in a system obtained by 
that transformation, the dynamics of the action variables I a and the radiation energy 
densities a a can be determined independently of the angle variables 9 a . In order to show 
this statement, we put several assumptions. 
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(i) The interaction energy of the scalar fields U' mt is analytic with respect to the dynam- 
ical variables I a , 9 a and is 2n periodic with respect to 9 a . Ui nt is bounded as 

^t-^l/l, (5.11) 
which implies that there exists a positive constant M such that 

|a 9/2 t/mt| < vM\I\, (5.12) 

where 

|/|:=53|J B |. (5.13) 

a 

(ii) As for the total energy density 

p = E^+Ei+^. ( 5 - 14 ) 

^— ' cr ^— ' a 4 

a a 

the action variables J a and the radiation energy densities cr a satisfy 

ci < a^mJa + < c 2 , (5.15) 

a a 

for some positive constants Ci, C2. 
Note that the evolution equations of the system can be written as 

a-fj-Ia = F a (I,a,9,a), (5.16) 
da 

a—a a = F a (I,a,9,a), (5.17) 
da 

a—9 a = -u a (I,a,a) + G a (I,a,e,a), (5.18) 
da e 

where F a F a G a are analytic with respect to the dynamical variables I a a a 9 a and 2ir 
periodic with respect to the angle variables 9 a . We say that the evolution equations is of 
the type 

(i) The averaged parts of F a F a G a are bounded as 

<F a > ~ a 2 \I\, (5.19) 

<F a > ~ a 3 \I\, (5.20) 

<G a > ~ a 2 , (5.21) 

where < A > implies the averaging over the angle variables 9 a : 

<a>: =WfI Vs * a (5 - 22) 

In case of the resonant case, it is prescribed that the averaging is performed with 
respect to the fast angle variables only [TOl]. 
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are 


bounded 


as 


F a 


~ e k a 2 \ 




F a 


~ e fc a 3 | 


I\, 


G a 


~ e k a 2 , 


i 



(5.23) 
(5.24) 
(5.25) 

where A implies the residual part after the averaging over the angle variables 9 a : 

A := A- < A> . (5.26) 

Under this notation, the following proposition holds. 

Proposition 1 Let k be some non-negative integer and consider the evolution equations 
as 

a ^/W = J?W(/W >(T W,flW,a), (5.27) 
a^af = ^(/W.aW.flW.a), (5.28) 
a^OW = i w f)(/(VW lfl ) + Gf»(/(V ( ^ (t) ,4 (5.29) 
Suppose that this set of evolution equations is of the type Ck, there exists a transformation 

/(*) = I^+u a k \l^ k+1 \^ k+1 \9^ +1 \a), (5.30) 
a (k) = (T (*+i) +u (*)(/(*+i) ><T (*+i) > 0(*+i) >a ) > (5.31) 
( fc ) = e ik ^ +v a k \l ik ^\a ik+1 \e ik+1 \a) ) (5.32) 

satisfying the following conditions: 

(i) wi^ v are analytic with respect to the dynamical variables li k+1 ^ aa +1 ^ a k+r> , 
are 2ir periodic with respect to the angle variables 6 a k+1 \ and are bounded as 

( fc ) ^ 



li 



(k) 



a 



|/ (fc+1) |, (5.33) 
e fc+1 a|/( fc+1 )|, (5.34) 
y( fc) ~ e fc+1 . (5.35) 

(ii) TTie evolution equations of the transformed variables &a +V> a k+1 ^ are of the 

type Ck+i and the changes of < F a > < F a > < G a > are bounded as 

A < F a > ~ e fc+1 a 2 |7 (fc+1) |, (5.36) 
A < F a > ~ e fc+1 a 3 |/ (fc+1) |, (5.37) 
A < G a > ~ e fe+ V, (5.38) 

where A < A > is defined by 

A < A >:=< A^ +1 \l^ k+1 \^ k+1 \9^ k+1 \a) >-< A^\l {k+1 \ a^ k+1 \ 6 {k+1 \ a) > . 

(5.39) 
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(For the proof, see Appendix A.) 

This proposition implies that we can make the part depending on the angle variables 
arbitrarily small by taking the original set of the evolution equations as the starting point 
k = and applying transformations given in the proposition iteratively. Therefore it 
can be expected that the evolution of the dynamical variables can be described by the 
truncated system obtained by discarding the angle variable dependent part with sufficiently 
good accuracy, if we take sufficiently large k. By estimating the errors produced by the 
truncation, we show that the above expectation is correct. We use the symbol A to 
represent the difference of a quantity for the exact system and a corresponding quantity 
for the truncated system. For A = (J , a a , 9 a ), the errors of the background variables are 
written as 

AA:= A- A tT , (5.40) 
and the errors of the perturbation variables are written as 

ASA := 5 A - 5A tr , (5.41) 

where A, 5 A represent quantities of the exact system and A tr , 5A tr represent quantities of 
the corresponding truncated system. For a function f(a), let us write 

f(a)=E(-a 2 ), (5.42) 

when f(a) is bounded as 

\f(a)\ < p(a) exp (-Aa 2 ) (5.43) 

for a polynomial of a: p(a), and for a positive number A. For a function /(a), let us define 
ll/(a)l|by 

||/(a)|| := sup |/(a')|. (5.44) 

l<a'<a 

When we write all the inequalities, it is prescribed that all the coefficients of order unity 
are omitted. The truncation error for the m-th order system can be estimated as follows. 

Proposition 2 A Let m be an integer larger than or equal to 2. For the m-th order system, 
the truncation errors of the background variables are given by 

\AI\ < E{-a 2 )e m , (5.45) 
\\Aa\\(a) < e m , (5.46) 
|A0| < a 2 e m -\ (5.47) 

and the truncation errors of the perturbation variables are given by 

\ASI\ < E(-a 2 )e m - 1 8A 1 (l), (5.48) 

\A5a\ < e^SA^l), (5.49) 

|A50| < a 2 e m - 2 5A 1 (l) + e" 1 " 1 exp (a 2 e m ) [a 2 e\59(l)\ + a 4 5A m (l)] , (5.50) 

where 

SA m (l) := + \5a(l)\ + e m \S9(l)\, (5.51) 



21 



under the initial conditions 



AJ(1) = A<j(l) = A9(l) = 0, (5.52) 
A<SJ(1) = A(J<t(1) = AS6{1) = 0. (5.53) 

(For the proof, see Appendix A.) 

By the transformation laws, the errors of the m-th order variables affect the original 
variables as shown by the next proposition. 

Proposition 2B The difference between obtained from A^ by the transformation 
laws and A[^ obtained from A^ by the same transformation laws has upper bound 



|A/(°)| < E(-a 2 )e m , (5.54) 

\AaW\ < e m , (5.55) 

|A0(°)| < a 2 e m -\ (5.56) 

and the corresponding difference as for the perturbation variables has upper bound 

\A5lM\ < S(-a 2 )e m - 1 ^ m) (l), (5.57) 

\A5a^\ < e^SA^il), (5.58) 



|A^ (0) | < a 2 e m - 2 54 m) (l) + e m - 1 exp(a 2 e m ) [a 2 t\56 {m \l)\ + a A 5A { ™\\)] , (5.59) 

under the same initial condition as in the previous proposition. 
(For the proof, see Appendix A.) 

According to the above proposition, we can conclude that we can make the trunca- 
tion errors as for the original variables arbitrarily small if we truncate the system at the 
arbitarily large m-th order system. 

From Proposition 1, we can see that the part independent of the angle variables of 
the evolution equations are shifted after the transformations reducing the part dependent 
on the angle variables. By the truncation, our system become much simpler than the 
original system. But it is still difficult to solve the truncated evolution equations with 
such correction terms because the evolution equations are complicatedly entangled with 
each other. In order to solve the evolution equations analytically, we want to discard such 
correction terms. The errors produced by discarding such corrections are evaluated in the 
following proposition. 



Proposition 3 The difference between the system with the correction terms produced by 
the transformation and the system obtained by discarding such correction terms is evaluated 
in the following way. As for the background variables, the discard errors are evaluated as 

\AI\ < E(-a 2 )e, (5.60) 
|Act| < e, (5.61) 
|A0| < a 2 , (5.62) 
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and as for the perturbation variables, the discard errors are evaluated by 



where 

under the initial conditions 



\A6I\ 


< 


E(-a 2 )e5B(l), 


\A5a\ 


< 


eSB(l), 


\AS6\ 


< 


a 2 5B(l), 


5B(1) 




\SI(l)\ + \Sa(l)\, 



(5.63) 
(5.64) 
(5.65) 

(5.66) 



AJ(1) = Acr(l) = A0(1) = 0, (5.67) 
A8I(1) = A8a(l) = AS6(1) = 0. (5.68) 

(For the proof, see Appendix A.) 

By the above propositions, it can be understood that we can obtain the information of 
the original system with sufficiently good accuracy by investigating the evolution equations 
by simply dropping the part dependent on the angle variables, because the errors produced 
by dropping are sufficiently small and in particular Act, A5a are bounded. The reason why 
these errors are mild is that the final state in which all the energy of the scalar fields is 
completely transferred into that of radiation fluids is the attracting equilibrium around 
which the perturbations do not grow. 



§6 Application of the LWL formula to the multicom- 
ponent reheating model 

In this section, we apply the LWL formula to the reheating where the energy of the multiple 
scalar fields is transferred into that of the multiple radiation fluids. The decay rate from 
the scalar field <p a to the radiation fluid p a is given by T aa . When the interactions between 
the scalar fields <p a - CAnt are neglected, the background quantities are solved as 

mj a = A a exp ^-jT a J^da-J—^ , (6.1) 

a a = B a + J^da-^^2T aa m a I a , (6.2) 

where A a , B a are integration constants. As long as we do not give the expression of 
the total energy density p in the integrals, the above solutions do not give any physical 
information of reheating. But it is difficult to solve the evolution equations in the form 
where the exact expression of p is explicitly described, because the evolution equations of 
I a , a a are complicated and highly nonlinear. Then we expect that the contribution to the 
integrations owes mainly to the period when the energy of the scalar fields is dominant, 
that is p can be approximated as 

A B . . 

" = ? + ?• < 6 - 3) 
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where 

A:=J2^a, B:=Y,B a . (6.4) 

a a 

We assume that the initial radiation energy density B is negligibly small B <C A. We 
assume that all Y a := J2 a ( a = 1; 2, • • •, N s ) are of the same order of magnitude. By 
substituting the above p expression to the solutions (I6.ip(l6.2l) . by expanding the solutions 
with respect to B around B = 0, we obtain 

mala = ^exp|-^a 3 / 2 } +7 r a ^aV 2e xp{-^a 3 / 2 ) + 0(5 2 ), 



4 4 \^ 

a p r = a 2_^Pa 



B I 1 + 2 G(2) " 2 G(1) | + (2 ) ° (5/3) W ^ i# + 0(52) ' (6 - 5) 



where 



r a :=^r aa , (6.6) 

a 

and G(t) is the Gamma function defined by 

POO 

G(t) : = / dxx^e'*, (6.7) 



which is convergent for t > 0. By taking the exactly homogeneous perturbation of (16.51) 
defined by 

D := [8A- 4j + 5b -tL) > ( 6 - 8 ) 



a B=0 



OA dB, 

we can obtain the Bardeen parameter in the final state a — > 00: 
I Dpi 



Ks) 2 ' 3 ^ I 1 + ^ G(2) 4 G «} / ^ (5/3) 



+ 4 1 3 a ^ r 2/3 ^ r 2/3 1 / ^ r 2/3 ' 



From now on, we name the expressions representing the final amplitude of the Bardeen 
parameter £jj n in terms of the initial perturbation amplitudes such as 5A a 5B a , S formula 
after S matrix in the quantum mechanics. As for the initial energy density perturbations 
of the scalar fields SA a , by adopting more physical parametrization introduced by 

ISA 6A a 6A b 

3T =: t> ab = ^ ~A^ =: ^° 6 ' ( 6 °) 

5A a can be written as 

6A a = 3A a ^ + J2^ L Va b . (6.11) 



A 

b 
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£ represents the adiabatic growing mode and r] ab represent the isocurvature modes. Since 
r] ab satisfy 



Vab 



-Vba 



Vab + Vbc = Va 



(6.12) 



the independent quantities are given by £ rji 2 rj 2 3 ■■■ Vn s -i n s - By using this parametrization, 
5A a dependent part of Cjj n is written as 



l J2 

8 L^i 



A a A b 



ab 



r 



2/3 



^2/3 



A 



-Vab, 



\- Ag 

2-^ r 2/3 ' 
a 1 a 



(6.13) 



where A D B implies that B is contained by A, that is A = B + • • -. From this S 
formula, we can conclude that for the adiabatic growing mode £, the Bardeen parameter 
is conserved, and that the initial entropy perturbations survive in the case that the decay 
rates are dependent on the scalar field <fi a from which the radiation energy comes, that 
is T a ^ Tj (a / i). In the case where multiple scalar fields exist, there is no reason 
why the perturbation has only the adiabatic component, and it is natural to think that in 
the perturbation the adiabatic components and the entropic components coexist. In such 
mixed cases, so called conservation of the Bardeen parameter does not hold and the above 
S formula gives useful tool for calculating the final Bardeen parameter. 

We can consider the case where the energy transfer rates T aa fluctuates [291]. which is 
called as the modulated reheating scenario [H]. For simplicity, we consider the one scalar 
field case. By taking the derivative of (16. 5ft with respect to T with B vanishing, we obtain 



C» D — - — 

s fin 6 r ' 



(6.14) 



which is well known formula derived in the paper [H], and where the coefficient is success- 
fully determined in this paper. 

We consider the influence of the resonant interaction between scalar fields on the final 
amplitude of the Bardeen parameter. We take the interaction U- m % into account by iteration. 
As the first order correction from the interaction term as 



we obtain 



where 



1 a 3 dU- mt d_ 
epV2 d6 a da a ' 



A (mJa) := exp |~7r o J^da-^—^ J^da^-exp S^jT a J^da-^ 
Since as the zeroth order approximation I a obeys (16.11) . we can write 



1 a 3 



e P 



1/2 



hit 



oe a de a 



a=l a 



3n/2 



exp < — 7r / da 



(6.15) 



(6.16) 



(6.17) 



(6.1* 



25 



where we assumed that U- mt contains an n-th order interaction term and that T is the 
appropriate sum of T a /2. By substituting the above expression and by expanding the 
correction term with respect to B around B = 0, we obtain 



i.m n 



mt 



a=l 



2 7 \ n_8/3 1 



7 




G(5/3,-n + 3,r a ,r)+(-^J ' f - ) (J/- r„)G'(r,/3.-n- 10/3.J-„.J-) 



1/3 



1 /2\ 2/3 /3\ 1/3 

- I - J G(5/3, -n + 7/3, r a , T) + I - J r a G(2, -n + 3, T a , T) 



1 (2 



2 V3 



2/3 



G(i,-n + 3,r a ,r) 



where 



G( ni ,n 2 ,T h T 2 ):= dxx n ^ exp (-T x x) / rf^- 1 exp {(r x - T 2 ) y} 



where 



7 2 



a 3/2 , x := 



7 2 



(6.19) 



(6.20) 



(6.21) 



41/23" ' 41/23 

The evaluation of the double Gamma function defined by (16.20p is treated in Appendix B. 
We consider the concrete example defined by 



U- int = A0 2 02, m 1 = m 2 . 
In this case, we use the independent variables defined by 



(6.22) 



0i 
h 



Po 



9 2 ~- 
-Pi 



Qo + qi 
h = Pi, 



(6.23) 



where (qo,Po) and (?i,Pi) are called fast and slow action-angle variables, respectively [TU1]. 
Because of the resonant relation satisfied by the masses of the scalar fields, the slow angle 
variable q± moves much more slowly than the fast angle variable go- The averaging over 
the slow angle variable q± cannot be justified, and therefore the slow action-angle variables 
(qi,Pi) can have evolutions. In the previous paper [HI], we investigated the influences 
of the resonant interaction on the evolution of the cosmological perturbations before the 
energy transfer from the scalar fields to the radiation fluids begins. According to this 
study, the slow action-angle variables can suffer from the instability near the hyperbolic 
fixed point in the phase space of the slow action-angle variables. Since the initial adiabatic 
perturbation £ and the initial isocurvature perturbation r\\ 2 are given by 



1 $Po 
3 Po 



and 



Pi5p - p Q 5pi 
(Po - Pl)Pl 



(6.24) 
(6.25) 
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respectively, the instability of the slow action variable pi has influence on the isocurvature 
mode. Therefore from (I6.13P in the case Fi ^ T 2 the instability of the action-angle variables 
survives in the final amplitude of the Bardeen parameter. Next we calculate the first order 
correction term (I6.19P in the present model (I6.22p . The present model has the hyperbolic 
fixed point at 

qi = l(2k + l), 2 Pl =p = c, (6.26) 

where k is an integer. At this hyperbolic fixed point, the first order correction to the final 
amplitude of the Bardeen parameter (16.191) is calculated as 

CL 3 T2A^-e C26qi(1) I ^ " ^ I / > ( 6 ' 27 ) 
where 

v : = (6.28) 
m o a o 

and the non-dimensional masses are scaled as mi = m 2 = 1. In the present model, as for 
the first order correction term also, in order that the slow action-angle variables instability 
has influence on the final Bardeen parameter, Y\ ^ T 2 is necessary. 

Until now, we evaluate I a a a by assuming that p is given by (16.31) . Now we evaluate 
the contribution to I a a a from the late stage of reheating when p is given by 




Plate 



'3\ 2/3 1 A 1 ' 3 v A 

a 1 a 

Such late stage of reheating begins at 

ai=d \2j G(5/3) 7 2/3 A 2/3 Z]^' ( 6 - 30 ) 

because at this ai, A/ a 3 is almost equal to pi a te- d is a numerical factor which can be 
assumed to be larger than unity. By using pi ate , ai, we can evaluate the contribution to p r 
from the late stage of reheating as 

i / q \ 2/3 I p 



Pr D G(5/3)^rfexp| -d 3 ' 2 G(h/?>f 



A 3/2 




X ^2/3/12/3 E r 2/3 ' (6.31) 
' b 1 b 

whose size is characterized by 

r l/e := dexp {-d 3/2 G(5/3) 3/2 } , (6.32) 

which is the ratio of the late contribution to the main early contribution to p r . The value 
of r// e is 0.43, 0.18 and 0.035 for d = 1,2 and 3, respectively. Therefore we can conclude 
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that the S formula which is derived by using (16.31) is rather good approximation to the real 
S formula. 

We consider the case where the decay rates T a depend on the radiation temperature T. 
In the high temperature limit T 3> m where m is the mass scale of the oscillatory scalar 
fields, the decay rate Y a depends on the radiation temperature [35' ]. When the decay 
product is the fermion, r a is given by 

T a = a a ^. (6.33) 

When the decay product is the boson, T a is given by 

T a = (3 a T. (6.34) 

According to the paper [35]] the reason is following. We consider the case where p r ~ T 4 
is sufficiently high. In the fermion case, the Pauli exclusion principle inhibits the decay of 
4> a into fermions since the fermions have already occupied the energy levels into which <p a 
would decay. In the boson case, the induced effect promotes the decay of 4> a into bosons 
since the bosons occupy the energy levels into which a decay. For simplicity, we consider 
the case where the radiation consists of one component. We interpret that the radiation 
temperature T appearing in (I6.33j) (16.341) is the temperature T a := [p r (a(r a ))] 1//4 at the time 
when the decay process proceeds given by 



a(r a ) := nln _ 2/3 - (6.35) 



A 1 / 3 

7 2/3 ra 

In the fermion case, by substituting T a defined above to (I6.33P it can be verified that 

10/9 



T,^3 = j^[Y,4k) ■ ( 6 - 36 ) 



Therefore we obtain 




(6.37) 
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E ( - -575 ) ^rW E ( 6 - 38 ) 



In the same way as in the fermion case, in the boson case we can obtain 

2/3 




r -1 



(6.39) 



1 ' 

ab 



The radiation temperature dependence of the decay rate T a affects how the isocurvature 
modes are transmitted into the final amplitude of the Bardeen parameter. 
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We consider the non-Gaussianity of perturbations. For simplicity, we assume that all 
the perturbations are generated by only one Gaussian variable, say locally homogeneous 
perturbed variable C(x) = C+SC(x) where SC(x) is spatially dependent Gaussian random 
variable. By assuming that p r (a, C(x)) oc C(x) a , from (I3.24p . the non-linearity parameters 
US] defined by 

C = Ci + pNLC! + 7^9nl(! + ■■■ (6.41) 

where (\ is the first order perturbation of the Bardeen parameter generated by one Gaussian 
random variable, are given by 

Inl = --, gsL = (6-42) 

In the reheating model with only one scalar field, the final radiation energy density is given 
by 

°V ~ ^73 pTF- ( 6 - 43 ) 

When the initial action variable A/m is the random Gaussian variable, by considering 
A(x) oc C(x) 2 from (|3.29p . non-linearity parameters are given by 

5 25 
Inl = -4, 9nl = — • (6-44) 

In the modulated reheating scenario, by assuming that the decay rate T(x) is proportional 
to <p(x)P, the non-linearity parameters are given by 

5 100 , N 

Inl = p 9nl = jjcp- (6.45) 

By observing the non-Gaussianity, we can determine whether the nonnegligiable Gaussian 
random variable lies in the action variable A(x)/m or the decay rate r(aj). 

§7 Application of the LWL formula to the multicom- 
ponent curvaton model 

In this section, we apply the LWL formula to the curvaton scenario where multiple weakly 
coupled massive scalar fields called curvatons decay into multiple radiation fluids some 
time later after the inflation has ended. In this curvaton scenario, the curvaton fields other 
than the inflaton fields driving the inflation are responsible for the origin of the cosmic 
structures. First, we assume that all Y a := J2 a ^aa (a = 1, 2, • • •, Ns) are of the same order 
of magnitude. 

First we consider the limit where in the initial time the curvaton fields energy A is small 
compared to the radiation fluids energy B: A <^ B. By substituting the p expression (16.31) 
to (16. ip (16.21) and by expanding it with respect to A around A = 0, we obtain 

mj a = A a exp |-^-^^a 2 | + 0(A 2 ), 

Dl/4 A 

a'p r = a*Y,P* = B + V2G(3/2)—J2-^ 2 +0(A 2 ). (7.1) 
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By taking the exactly homogeneous perturbation, that is D operation around A = 0, we 
obtain 

A ^ 1 Dffi _ 1 SB y/2 1 v SA a 

Cfi„ - ~ A — - ~ ~ A — + — G(d/2J^ 1/2/?3/4 2^ ^172- (7.2) 



4 p 4 5 



7 l/2 £3/4 ^ pi/ 



Next we consider the case where the energy densities of the curvatons are large compared 
with those of radiation fluids when the energy transfer from the curvatons to the radiation 
fluids proceeds. In this case, the exponent of (16.11) is written as 



7r a I da— 

where 

As for x(T a ) defined by 



i 7r a (2 



*-x^-2Bx^ + -B^ 
x p l ' z a A A \6 3 



we obtain 



x(T a ) = x (T a ) 



7^ 



1 + 2 



x := Aa + B. 

{x{Y a )-B}/A 

da 



p l l 2 a 



B 



x (T a ) 3 V^o(ra) 



B 



3/2 



+ 



B 2 



^o(r a ) 



where 



/3\ 2/3 /A 2 \ 2/3 



(7.3) 

(7.4) 
(7.5) 

(7.6) 
(7.7) 



The expansion parameter B/xo(T a ) implies the ratio of the energy of radiations to the en- 
ergy of the curvaton a when the energy transfer proceeds. By using x(T a ), we approximate 
the exponential function by the step function: 



e ^{- lTa l da pk} 



(x(T a ) - x) 



(7- 



By using this approximation, we obtain the radiation energy density in the a — > 00 limit: 



Pr 



V- 1 D 2 /3\ 5/3 ^ 1/3 v- 



A n 



ry2/3 p2/3 

a 1 a 



B 



(7.9) 



By taking the exactly homogeneous perturbation, we obtain the S formula: 

1 5 A* 1 5B B 



fill 



+ 



4 A* 2 A, 2A\ 



<L4* + ■ 



where 



A * : =5 2 



2 /3\ 5/3 A 1 / 3 ^ A 



t^2/3' 



(7.10) 



(7.11) 
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We define more physical parametrization by 



8A a = 3A a ^ + 3A aV + J2^^Vab, (7.12) 

b 

5B a = 4iu+E^ir^- ( 7 - 13 ) 

In particular, since 

5 A = 3Af + 3Ar], 5B = 4B£, (7.14) 

£ implies the adiabatic growing mode and i] means the isocurvature mode between the 
total curvatons and the total radiations. By using this parametrization, the S formula can 
be rewritten as 



A 



+ 



[^-^)^r^/J2^k + ---- ( 7 - 15 ) 

* 7 ab \ La L b / a 1 a 



In the most simple curvaton scenario of one curvaton field and one radiation fluid, the 
empirical S formula was obtained from the numerical calculation [7]]: 

Cfin = r(p)r], (7.16) 
0.924 x " J - 21 



r(p) = 1 -[ 1 + T2A P I (7 - 17) 
where 

P := 7 l/ 2r i/ 2E 3/4 - ( 7 - 18 ) 

In the limits p <C 1, p ^> 1, this empirical S formula gives 

Cfi„ = 0.924^ p < 1, (7.19) 
/ 1.44 \ 

Cfin = ( 1 - ^2lJ ^ P > 1, (7-20) 

respectively. Our analytical results (17. 2ft ( 17. 15ft give 

Cj n D 0.940p77 p < 1, (7.21) 

h ( 2.54\ 

d D (l - J^) V V » 1, (7.22) 

respectively. In the case p < 1, the empirical formula and our analytic result agree with 
good accuracy In the case p ^> 1, our analytic result is obtained by rather rough treatment 

approximating the exponential function by the step function. But our analytic S formula 
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agree well with the empirical formula. For reference, for p ^> 1, according to our method, 
a more precice calculation gives 



(I 3 



10 
~3~ 



2/3 



10 1 10 (2 



V 



4/3 + 3 p 2 + g 



4/3 



P' 



8/3 



1 - 2.54^L + 3.33i_ + 0.647 1 



P 



1.33 



V 



p 



2.67 



P' 



3.33 



o 



10/3 



P 



(7.23) 



The errors between the above formula and the empirical formula are 1.3 percent, 1.1 percent 
for p — 10, p — 5, respectively. 

In section 6, until now in section 7, we assumed that all decay rates r a are of the same 
order of magnitude. As application example of the formulae (17.11) (17.91) . we consider the 
reheating where the decay rate of <fri is much larger than the decay rate 02; Ti 3> 1^2- Just 
after the scalar field 0i decays, by using the result of section 6, we obtain 



4 

a p r 



A 1 / 3 A x 

7 2/3 p2/3 • 



a 3 p s 



An. 



(7.24) 



By regarding a 4 p r in the above as B, we can use the formulae (17.11) (17.91) of the curvaton 
scenario. In the case A 2 /A x (I^/Ti) 1 / 2 , by using ( 17.11) . the final radiation energy density 
is calculated as 

A 1 '* A l 1 A 1 / 12 A\ /A A 2 , x 

1 " (7.25) 



4 

a p r 



+ 



^2/3 p2/3 1 ^2/3 pl/6pl/2 ' 

In the case A^ 3 '/ A 1 / 3 A x ^> (I^/Tx) 2 / 3 , by using (17.91) . the final radiation energy density is 



calculated as 



4 

a p r 



A 1 / 3 A x 1 At /3 



^2/3 p2/3 1 ^2/3 p2/3 ' (7.26) 

By substituting the above two expressions to (13.241) and by expanding with respect to 
Gaussian random perturbations, we can obtain the Bardeen parameter of arbitrary order. 



§8 Discussion 

In this paper, we constructed the LWL formula expressing the long wavelength limit of evo- 
lution of cosmological perturbations in terms of the corresponding exactly homogeneous 
perturbations in the most general scalar-fluid composite system. We determined the correc- 
tion term which corrects the difference between the long wavelength limit of cosmological 
perturbations and the exactly homogeneous perturbations, and we showed that the cor- 
rection term contributes the well known adiabatic decaying mode. It was pointed out that 
when we extract the long wavelength limits of evolutions of cosmological perturbations from 
the exactly homogeneous variables, the use of the scale factor a as the evolution parameter 
is more useful. The scalar-fluid composite system whose LWL formula is constructed in 
this paper can be used to discribe the early stage of the universe such as reheating after 
inflation and the curvaton decay in the curvaton scenario, when the fluid is assumed to be 
radiation. In this paper, the LWL formula is applied to the most general case of reheating 
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and of the curvaton decay containing the multiple scalar fields and the multiple radiation 
fluids, and the S formulae representing the final amplitude of the Bardeen parameter in 
terms of the initial adiabatic and entropic perturbations are constructed. In case where for 
different a, the value of the decay rate T a is different; that is T a ^ Tb for a ^ b, the initial 
isocurvature modes survive in the final amplitude of the Bardeen parameter. 

We discuss the non-linear generalization of the LWL formalism. Recently the gradient 
expansion has been discussed as the method for investigating the evolutions of non-linear 
perturbations on superhorizon scales [27] ] . In the lowest order of the gradient expansion, 
in the zero curvature slice da(t,x)/dx l = 0, the evolution equation of the scalar quan- 
tity T has exactly the same form as that of the exactly homogeneous equation of T [25] 
]. But the coefficients of the evolution equation are spatially dependent, therefore this 
evolution equations describes the non-linear superhorizon scale inhomogeneities. By using 
the solution of the exactly homogeneous system with the scale factor as the evolution pa- 
rameter T(a,C), the solution of the locally homogeneous evolution equation is given by 
T(a,C(x)), where C(x) = C + 6C(x) is spatially dependent solution constant. The spa- 
tially dependent perturbation part of T(a, C(x)) is given by Taylor expanding with respect 
to C(x) = C + 5C(x); 

°° 1 d k T 

T(a, C(x)) = T(a, C) + £ E ' ' ' E k\ dC ■ ■ ■ dC 5 ° ai{x) ' ' ' 6Cak{x) > (8 ' 1} 

k=l ai a k ' ai ak 

whose first order perturbation part agrees with our linear perturbation formula with ne- 
glecting the adiabatic decaying mode; DT^ = (dT/dC) a . Therefore T(a,C(x)) is the 
non-linear generalization of our linear perturbation variable DT. When P r = p r /3 as in 
the final state of reheating or the curvaton decay, T = lnp r /4 is the non-linear generaliza- 
tion of the Bardeen parameter (. Since in this paper we determine p r in the final state of 
reheating or the curvaton decay with the scale factor as the evolution parameter, we can 
obtain the information of the non-linear evolution and the non-Gaussianity of perturbations 
which fluctuate spatially on superhorizon scales. 

Our evolution equations have arbitrary functions S which describe the energy transfer 
between scalar fields and perfect fluids. The source functions S can be determined from 
the microscopic dynamics between the coherently oscillating scalar fields and radiation, 
concretely speaking, by path integrating out the fields constituting the radiation and inter- 
acting with the coherently oscillating scalar fields in the effective in-in action [33]- When 
the scalar fields oscillate coherently, the source term such as S = Tcj) is important in order 
that the energy is transferred from the scalar fields into radiation effectively. As shown in 
the paper [35]], Y is given as the function T = T(<p, 0, p r ) where p r is the energy density of 
radiation. As seen in the section 6, when T is a function of the scalar quantities which are 
closely related to reheating process, the functional form of T affects how the initial isocur- 
vature components are converted into the adiabatic component such as the final amplitude 
of the Bardeen parameter, but it does not affect the evolution of the initial adiabatic grow- 
ing mode. On the other hand, the scenario where the energy transfer is controlled by the 
scalar quantities not related to reheating is considered as the modulated reheating scenario 
[2"51 ] [H]. As such scalar quantities, we can choose flat direction scalar field which does 
not govern the energy of the universe but fluctuates of the order of the Hubble parameter 
during the inflationary expansion, or the scalar field written in terms of the slow action 
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variable which suffers from the hyperbolic instability due to the resonance of the masses 
of the scalar fields during the oscillatory stage PHI]- As seen from the formula (16.141) the 
fluctuations of such modulating scalar quantities are imprinted on the p r fluctuation. 

We consider the system where multiple oscillatory scalar fields and multiple radiation 
fluids interact. As for the system without radiation, its evolution of cosmo logical pertur- 
bations have been investigated in detail [9]], [10]]. The system of multiple scalar fields only 
can be written in terms of Hamiltonian form. Any Hamiltonian system obeys the Liouville 
theorem, that is, volume occupied by group of orbits are invariant. Therefore according to 
the LWL formula, in the stable case the perturbations do not grow and in the unstable case 
the same numbers of growing modes and decaying modes appear because of the squeezing 
of the phase space volume. The former case occurs in the case where masses of scalar 
fields are incommensurable and near the elliptic fixed points in the case where masses of 
scalar fields are commensurable. The latter case occurs near the hyperbolic fixed points 
in the case where masses of scalar fields are commensurable. Then we include dissipative 
interaction with radiation. In this case our system is not a Hamiltonian system and it does 
not obey the Liouville theorem. In this dissipative system, in addition to two possibilities 
mentioned above we can expect the third possibility where group of orbits are attracted 
into the attracting set. The final state where all the energy of the scalar fields is transferred 
into that of radiation fluids is the attracting equilibrium. Around the attracting set, the 
adjacent orbits come nearer and nearer, therefore the LWL formula tells us that all the 
perturbation modes are stable, that is, converge into some constants or decay. It is useful 
to investigate how the behavior of cosmological perturbations around the hyperbolic fixed 
points is changed due to the dissipative interaction with radiation under the spirit of the 
LWL formula. In this line of researches, we will have new explanation of the backreaction 
which suppresses the instability due to the resonance. In the future publication, we will 
return to this problem. 
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§A Proofs of the propositions in §5 
A.l Technical Lemmas 

Technical Lemma 1 Under the assumption Ii5.15\) . for \f\ < 1, 



and for \ f\ < 



df 


< a, 


df 
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< -, 


dl 




da 


a 


df 


<1, 


df 


<- l J l 


dl 




da 


— a 



(A.l) 



(A.2) 
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Technical Lemma 2 For the general physical quantity A(I, a, 9) with Fourier decompo- 
sition as 

A = Ao(I, o-) + J2 MI, ex P ( ik ■ 0). ( A -3) 

the solution to the first order partial differential equation 

d 

^a^S = A-<A>, (A.4) 

is given by 

S = {A}, (A.5) 

where 

w^E^ «•><*■«)■ ( A6 > 

A. 2 Proof of Proposition 1 

In order to make the notation simple, we omit the superscript (k) and replace (k + 1) with 
(1). By substituting the transformation laws of I a 9 to the evolution equations of I a a a 
6 a , we obtain 

<F a > +F a (I, a, 9, a) - F a (/«, a™ , 9^, a) - 

oa 

- p(i) ■ dUa F w ■ dUa F w ■ dUa (A 7^ 

- Fa + W ] b d$ " ' ( } 

< F a > +F a (I, a, 0, a) - F a (J< VW, a) - 



(9a 



du„ „m <9w„ „ m du r 



p(l) , ^ p(l) , p(l) , ""a r (l) ( . R x 

F « +^7(T)^ ' (A ' 8) 



dr b L > daf w\ 



and 



1 (t \ 1 ^r(i) (i) \ 1 duJa 1 <9^ a 
-w (J, <r, a) uj a {F >,a K >,a) m u h j^up 

e e e dl b ' e da\ ' 

+ <G a > +G a (I, a, 9, a) - G a (J« cr (1) , 0« a) - 

oa 



= ^ + ^< + ^> + ^Gl», (A.9) 
when we choose u a u a v a as 



dl^ b do-V * 89^ 



u a = e{F a } ~ e k+1 \I\, (A.10) 

u a = e{F Q }~e fe+1 a|/|, (A.ll) 

v a = l^^+^up + eGal-^. (A.12) 
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As for AF a AF a AG a defined by 



AF a 
AG„. 



F^ ] - < F a >, 
F«- < F a >, 
G«-<G a >, 



applying the mean value theorem to ( 1A.7I) ( 1A.8I) ( 1A.9h gives 



AF a + e k+1 AF b + e K+1 \I\AF + e« +x \I\AG 



k+l\ 



AF a + e k+L aAF b + e + a\I\AFp + e k+l a\I\AG b = 
AG a + e k+1 aAF b + e k+1 AF p + e k+1 AG b = e fc+ V, 



a z \I\ 



a 6 I 



(A. 13) 
(A. 14) 
(A.15) 



(A.16) 
(A.17) 
(A.18) 



AF a 


= e k+l a 2 


AF a 


= e k+1 a 3 


AG a 


= e k+1 a 2 



where all coefficients of order unity are omitted. By solving the above three equations, we 
obtain 

(A.19) 
(A.20) 
(A.21) 

AF a , AF a , AG a are decomposed into the angle variables independent parts A < F a >, 
A < F a >, A < G a > and the angle variables dependent parts Fa F^ Ga- 

A. 3 Lemmas and the preparatory propositions 

Lemma 1 The solution to the differential equation 



a^-A = -Xa 2 A + E{-a 2 )B{a) 
da 



where X is a positive constant, is bounded as 

A 



|A(a)|<exp(--a 2 )|A(l)| + J E;(-a 2 )||i?||(a) 



(A.22) 



(A.23) 



Lemma 2 When B satisfies 

d 



da 



-B 



< E{-a 2 )\B\ +E(-a 2 ){C+ \\B\\(a)} 



where C is a positive constant, for an arbitrary a > 1 

ai 



B a) < 



(C +11511(0!)) 



(A.24) 



(A.25) 



ai — 1 

where a\ is a constant satisfying ai > 1. For example, by putting a x = 2 we obtain 

\\B\\(a) <C +\\B\\(2). (A.26) 
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Proof By solving the differential equation, we obtain 



\B(a)\ < 15(1)1 + J daE(-a 2 ){C + \\B\\(a)}. 



Since for a satisfying 1 < a < oi, 



and for a > a% 



\B(a)\<C+\\B\\( ai ), 



\B{a)\<C+\\B\\{ ai ) + E{-$)\\B\\{a) ; 
then for an arbitrary a > 1 

||S(a)|| < H^IKa!) + £7(-a?)ll^ll(«) 



(A.27) 

(A.28) 
(A.29) 

(A.30) 



whose right hand side is an increasing function of a. Since E(—a\) < 1/fli, we obtain 



Lemma 3 When 81, 5a satisfy 

da 

-8a 
da 



< a\8I\ + \5a\ + \A\, 

< a 2 \8I\ +a\8a\ + \B\ 



the following inequality hold: 

||HI(2)<i5/(i)| + |Mi)i + HI(2) + P||(2). 

Proof We consider the differential equation as for a\8I\ + \8a\: 

4- (a\8I\ + \8a\) < a (a\8I\ + \8a\) + a\A\ + \B\. 
da 

Then we obtain 



(A.31) 
(A.32) 

(A.33) 
(A.34) 



\\8a\\(a) < exp ^-a 2 ){\5I(l)\ + \5a(l)\ + a 2 \\A\\(a) + o||S||(o)}. (A.35) 
We put a = 2. 

Proposition Apl For the m-th order system, for the background quantities, the following 
inequalities hold: 

|/| < exp(-a 2 ), \a\ < 1. (A.36) 



Proof We solve the evolution equations given by 

< -a 2 \I\+a 2 e m \ 



d T 

a—1 

da 

d 

a— a 

da 



< a 6 \I\ + a 6 e 



Zjm I 



(A.37) 
(A.38) 
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Proposition Ap2 Let m be an integer larger than or equal to 2. For the m-th order 
system, for the perturbation quantities, the following inequalities hold: 

(1), (A.39) 

(A.40) 



\5I(a)\ 


< 


5a\\(a) 


< 


\69(a)\ 


< 



1 m ^ - 1 - ) i 

.2 



,2,m\ 



a 



(A.41) 



where 8A m {l) is defined by 115. 51}) . 

Proof The perturbation variables satisfy the evolution equations as 

4-51 = -aSI + a\I\5o- + ae m \I\56, (A.42) 
da 

4 L 5a = a 2 5I + a 2 \I\5o- + a 2 6 m \I\56, (A.43) 
da 

A 2 

—56 = —5I+-5a + ae m S9, (A.44) 
da e e 

where |/| means a function bounded by M\I\ for a positive constant M. It is important to 
notice the coefficient of 51 in ( 1A.42I) is negative. Although in (1A.42I) . terms such as a 2 \I\5I 
is contained, such terms can be neglected because 

daa 2 \I\<\. (A.45) 



a 3 ,, „„, , , a 2 



By substituting the estimation obtained from ( 1A.44I) 

\59\ < exp (a 2 e m ) ||^(1)| + j\\SI\\(a) + ^\\Sa\\(a) j , (A.46) 

into ( 1A.42I) . by applying Lemma 1, we obtain 

\6I(a)\ < E(-a 2 ) {\5I(1)\ + e m |^(l)| + \\5a\\(a)} . (A.47) 
By using the above inequalities in (1A.43I) and applying Lemma 2, we obtain 

||<5a||(a)<| ( 5/(l)|+ e m |^(l)| + ||5a||(2). (A.48) 
On the other hand, by applying Lemma 3, we obtain 

\\6a\\{2)<5A m (l). (A.49) 

Then we can prove the results. 
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A. 4 Proof of Proposition 2A 

The evolution equations of AA where A = (J, a, 9) are given by 

= -aAI + a\I\Aa + ae m \I\, (A.50) 

da 

-^-Aa = a 2 AI + a 2 \I\Aa + a 2 e m \I\, (A.51) 
da 

A 2 

--A9 = —AI + -Aa + ae m . (A.52) 
da e e 

In the same way as the proof of the Proposition Ap2, we obtain 

|A/| < £(-a 2 )(|A/(l)| + |Acx(l)|+e m ), (A.53) 

\\Aa\\(a) < |A/(l)| + |Aa(l)|+e"\ (A.54) 

\A9\ < |A^(1)| + ^(|A/(1)| + \Aa(l)\ + e m ) , (A.55) 

In the above estimations, we put AA(1) = 0. 

Next we consider the evolutions of ASA where A = (I, a, 9). We take differences 
between the upper equations and the lower equations. As for A5I 

^-51 = -aSI + a\I\5a + ae m (51 + \I\5a + \I\59) , (A.56) 
da 



as for A5a 



and as for A59 



d ST 

da 


-a5I tT + a\I t]: \5<j tT , 




d 
da 


a 2 5I + a 2 \I\5a + a 2 e m 


(5I+\I\5a+\I\59), 


d A 

da 


a 2 5I tr + a 2 J tr \5a tI , 




da 


= ^ § I+-Sa + ae m 

£ € 


(a5I + 5a + 59) , 


^69 tT 
da 


a 2 a 
= — dl tr + -dcr tr . 
e e 





(A.57) 

(A.58) 
(A.59) 

(A.60) 
(A.61) 

By taking into account the fact that in the first two terms on the right hand sides the 
coefficients depending on /, a, not on 9 are multiplied, and using the estimations of 5A, 
AA where A = (I, a), we obtain 

5I-5I tv = ASI + 5I(aAI + Aa) 

= A5I + E(-a 2 )5A m (l)e m , (A.62) 

1 5a - I tr 5a tI = I tT A5a + 5a(AI + |/ tr |Afr) 

= E{-a 2 )A5a + E{-a 2 )5A rn {l)e m , (A.63) 

5a - 5a tT = A5a + 5a(aAI + Aa) 

= A5a + 5A m (l)e m . (A.64) 
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By using estimations of 5 A where A = (I, a, 9), we obtain 

1 



51 + \I\8a + \I\80 = -£?(-a 2 )5Ai(l), (A.65) 



e 



ao7 + 5a + 56 = exp (a 2 e m ) + ^<L4 m (l)) . (A.66) 



Therefore we get 



-^-A5I = -aA5I + E(-a 2 )A5a + e m - 1 E(-a 2 )(5A 1 (l), (A.67) 
aa 

-f-Aoa = a 2 A5I + E(-a 2 )A5a + e m - 1 E(-a 2 )5A 1 (l), (A.68) 
eta 

d A m fl2 a rr ° A r 

— A5# = — AS I + -A5a 

da e e 

+ae m ~ 1 exp (a 2 e m ) (e|<?0(l)| + a 2 <L4 m (l)) . (A.69) 

In the same way as the proof of the Proposition Ap2, we obtain 

\A5I\ < E{-a 2 ) {\A5I{1)\ + \A5a{\)\ + e m - l 5A 1 {l)) , (A.70) 

\\A5a\\(a) < \A5I(1)\ + \A5a(l)\ + e^SA^l), (A.71) 

|A<50| < |A(J0(l)| + ^(|A(J/(l)| + |A(J(7(l)| + e m - 1 (JA 1 (l)) 

+e m " 1 exp (a 2 e m ) (a 2 e|o"0(l)| + a A 5A m (l)) . (A.72) 

In the above estimations, we put A<L4(1) = where A = (I, a, 9). Then we complete the 
proof. 

A. 5 Proof of Proposition 2B 

From Proposition 1, we obtain 

j(o) = jM + e |/M| ; (A.73) 

a (0) = a (m) +ea \ I (m)^ ( A>74 ) 

0(0) = (m) +e> ( A>75 ) 

where |/( m )| means the function of A^ where A = (I, a, 9) bounded by M|/( m )| for a 
positive constant M. As for AA where A = (7, a, 6>), we obtain 

|A/ {0) - A/ (m) | < e (|A/ (m) | + |/ {m) ||A(x (m) | + |/ (m) ||A0 (m) |) 

< e m E{-a 2 ), (A.76) 
|Aa {0) - A(j (m) | < ea(|A/ (m) | + |/ {m) ||Aa (m) | + |/ (m) ||A^ (m) |) 

< e m £(-a 2 ), (A.77) 
| A0(°) - A0 (m) | < e (a| A/ (m) | + | Ao- (m) | + | A0 (m) |) 

< e m a 2 . (A.78) 
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By using the estimations of Proposition 2A and the above evaluations, we obtain the results 
of the former part. 

Next we consider A5A where A = (I, a, 9). By taking the variations of the transforma- 
tion laws, we obtain 



5 /(0) _ S j(m) = £ ^j(m) + |/M | 5(J (m) + | 7 (m) ^(m)) ^ 

5a (0) - 5(j (m) = ea (5/ (m) + |/( m )|5 a (m) + \I^\56 {m) ) , 

50 (O) _ S9 (m) = e ( a(J/ (m) + 5(J (m) + ^(m)J ? 



(A.79) 
(A.80) 
(A.81) 



where the coefficients are the functions of A^ where A = (I, a, 9). We take the differences 
of the transformation laws of the exact variables 5 A and those of the truncated variables 
5A tI . By using 





- 


= A5I {m) + 


5/ (m) (aA/ (m) + Aa (m) + A9 (m) ), 


(A.82) 


I5a im) - 




= 7 tr A5(j (m) 


+ 5a (m) (A/ (m) + / tr Aff W + / tr A# (m) ), 


(A.83) 


5cr {m) 


£ ( m ) 


= A&r (m) + 


5a {m \aAI {m) + Aa (m) + A# (m) ), 


(A.84) 


5 g(m) 


- ^ 


= A59 {m) + 


5# (m) (aA/ (m) + A(T (m) + A9 (m) ), 


(A.85) 


159 {m) - 




= I tI A59 (m) 


+ ^( m )(A/ (m) + J tr Aa (m) + J tr A# (m) ), 


(A.86) 



and by using the estimations of AA^ and 5A^ where A = (J, a, 9), we obtain 

A5I {0) - A5I im) = e m - 1 E(-a 2 )5A 1 (l), (A.87) 
ASa® - ASa™ = e m ' 1 E(-a 2 )5A 1 (l), (A.88) 
A59^-A59 {m) = exp (a 2 e m )e m " 1 (a 2 e\59(l)\ + a 4 5A m (l)) . (A.89) 

By using the estimations of Proposition 2A and the above evaluations, we obtain the results 
of the latter part. We complete the proof. 



A. 6 Proof of Proposition 3 

As for A A where A = (I, a, 9), we can obtain the results by putting m — 1 in the proof of 
Proposition 2A. 

Next we consider the evolutions of A5A where A = (I, a, 9). We take differences 
between the upper equations and the lower equations. As for ASI 



d_ 

da 



81* 



d ST 

-1-0 1 Q 

da 



as for A5cr 



5a, 



tr 



d_ 

da 

d A 

da 



-a5I tl + a|/ tr |5<7 tr + ae (5I tl + |itr|<fo"tr) , 
—aSI n + a\I n \5a n , 

a 2 5I tT + a 2 |/tr|<W + a 2 e (5I tI + |itr|<W) , 
a 2 5I n + a 2 \I n \5cr n , 



(A.90) 
(A.91) 

(A.92) 
(A.93) 
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and as for A56 

d a 2 a 

—59 tr = —5I tr + -5a tr + ae(a5I tr + 5a tT ), (A.94) 
da e e 

d 59n = ^ 5In+ a 5(Jn (Ag5) 
da e e 

In the above equations, the subscript tr implies that in the present system, the angle 
variables dependent parts which have been made sufficiently small by the transformations 
defined in the proof of Proposition 1 have already been truncated, and the subscript n 
means the further neglection of e-order corrections produced by such transformations. We 
obtain the evolution equations: 

4- AS I = -aA5I + E(-a 2 )(A5a + e5B(l)), (A.96) 
da 

-^-A5a = a 2 A5I + E(-a 2 ) (ASa + eSB(l)) , (A.97) 

— AbQ = —A5I+-A5a + a5B(l). (A.98) 
da e e 

In the same way as the proof of the Proposition Ap2, we obtain 

\A6I\ < E{-a 2 ) {\A8I{1)\ + \A8a{\)\ +eSB(l)), (A.99) 

||A<H(a) < |A5/(l)| + |A5<r(l)|+e5£(l), (A.100) 

a 2 

\A56\ < \A56(l)\ + — (\A5I(l)\ + \A5a(l)\+e5B(l)). (A.101) 

By putting ASA(1) = where A = (J, a, 9) in the above inequalities, we obtain the results 
of the latter part. We complete the proof. 



§B Evaluation of the gamma-like function 

In the present paper, we often have to evaluate the integrals defined by 



G(t,T) := / dxx l - x e- v \ (B.l) 

J x 

where x is defined by 

2 7 

X ° := 3AV^' (R2) 
where 7 is assumed to be sufficiently small. By expanding with respect to the small 
parameter x by the partial integration, we obtain the evaluations as follows. For t > 0, 

G(f,r) = lG(f) = 0(l), (B.3) 

for t = 0, 

G(0,T) = -mx + O(l), (B.4) 
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for t < and t is an integer, 

G(t, T) = -^4 + 0(4 +1 , Inzo), (B.5) 

and for t < 1 and t is not an integer, 

G(t,T) = - 1 1 x t + O(x t +\l). (B.6) 

G(t) is the well known Gamma function. 

Next by using the above evaluations, we evaluate the integrals defined by (16.201) which 
appear when we evaluate the effects of the interactions between scalar fields on the final 
radiation energy density p a = a a /a 4 . By expanding with respect to the small parameter 
Xq by the partial integration, we obtain the evaluations as follows. We assume that n\ > 0. 
For n 2 > 0, 

G(n 1 ,n 2 ,r 1 ,r 2 ) = 0(l), (B.7) 

for n 2 = 0, 

G{m, o, r 1; r 2 ) = --^ \nx G( ni ) + 0(1), (b.s) 
1 1 

and for n 2 < 0, 

G(m, n 2 , F 1} T 2 ) = -l^x^GK) + 0« + ™ 2 , x™ 2+1 , lnx ). (B.9) 

n 2 1 1 

Finally we evaluate the incomplete Gamma function defined by 

poo 

G(t; Xl ) := / dxx^e-*, (B.10) 

for large xi. By partial integration, we obtain 

G(t; xi) = x\~ l e- Xl + 0(a;'f V X1 ), (B.ll) 

for sufficiently large x%. 
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